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We present a new set of three-body interaction models based on the Bruch-McGee (BM) potential 
that are suitable for the study of the energy, structural and elastic properties of solid "^He at high 
pressure. Our ah initio three-body potentials are obtained from the fit to total energies and atomic 
forces computed with the van der Waals density functional theory method due to Grimme, and 
represent an improvement with respect to previously reported three-body interaction models. In 
particular, we show that some of the introduced BM parametrizations reproduce closely the exper¬ 
imental equation of state and bulk modulus of solid helium up to a pressure of ~ 60 GPa, when 
used in combination with standard pairwise interaction models in diffusion Monte Carlo simula¬ 
tions. Importantly, we find that recent predictions reporting a surprisingly small variation of the 
kinetic energy and Lindeman ratio on quantum crystals under increasing pressure are likely to be 
artifacts produced by the use of incomplete interaction models. Also, we show that the experimental 
variation of the shear modulus, C 44 , at P < 25 GPa can be quantitatively described with the new 
set of three-body BM potentials. At higher pressures, however, the agreement between our Ca 
results and experiments deteriorates and thus we argue that higher order many-body terms in the 
expansion of the atomic interactions probably are necessary in order to better describe elasticity in 
very dense solid ^He. 
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I. INTRODUCTION 

The electronic structure of a single ^He atom is among 
the simplest in the periodic table of elements. Likewise, 
the atomic interactions in liquid and solid helium can 
be reproduced accurately with simple analytical func¬ 
tions that solely depend on the distance between par¬ 
ticles taken in pairs. Examples of successful "^He-'*^He 
interaction models include the Lennard-Jones and Aziz- 
type semiempirical potentialsYet, under conditions 
of large pressures and strain deformations the interparti¬ 
cle interactions become more complex due to the strong 
electronic repulsion experienced by neighboring atoms. 
Consequently, pairwise potentials, which work reason¬ 
ably well under near-equilibrium conditions, turn out to 
be unreliable. This is, for instance, the case of the Aziz-II 
potential)^ which at high pressure provides too repulsive 
atomic forces and a significant overestimation of the ^He 
molar volume and bulk modulus.— 

A recently proposed straightforward way to correct for 
such modeling drawbacks consists in modifying the repul¬ 
sive part of standard pairwise potentials by means of an 
exponential attenuation factor.— This possibility has al¬ 
ready been explored in highly compressed solid and 
molecular hydrogen^ with quantum Monte Carlo simula¬ 
tions, producing equations of state which are in very good 
agreement with experiments. Nevertheless, the use of 
modified pairwise potentials in very dense crystals poses 
a series of issues and open questions. For instance, a 
surprisingly small variation of the kinetic energy upon in¬ 
creasing pressure have been reported in works @ and Q , 
and, owing to the lack of experimental data in the ther¬ 


modynamic regime of interest, it remains to be demon¬ 
strated whether such predictions can be fully ascribed to 
genuine quantum nuclear effects or not. Also, pairwise 
potentials are in general not recommended for the study 
of elasticity in hep crystals at high pressure since they 
inevitably lead to null values of the Cauchy discrepancy 
(defined as the difference between the two elastic con¬ 
stants C12 and (744), in contrast to what is observed in 
experimentsi^^— 

An alternative route to improve the description of 
quantum solids under extreme stress-strain conditions 
is to consider higher order terms, beyond pairwise ad¬ 
ditivity, in the approximation to the atomic interac¬ 
tions. In this context, several three-body interatomic 
models have already been proposed like, for instance, the 
Axilrod-Teller (AT), Bruch-McGee (BM), and Cohen- 
Murrel (CM) potentials— However, improvements 
resulting from the use of those three-body interaction 
models so far have been reported to be only marginal. 
For instance, three decades ago Loubeyre claimed, based 
on the outcomes of self-consistent phonon and classical 
Monte Carlo simulations, that the BM three-body inter¬ 
action could bring into good agreement calculations and 
experiments performed on the equation of state of solid 
helium up to ~ 60 GPa— i However, Chang et alr^ have 
shown more recently that when either the BM or CM 
three-body potentials are considered in quantum Monte 
Carlo simulations the resulting ^He molar volumes are 
significantly underestimated, already at few GPa. Simi¬ 
lar discouraging results have been reported also by other 
authors who have employed analogous three-body inter¬ 
action models— 
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In this article, we present new work done on the model¬ 
ing of three-body interactions in highly compressed solid 
helium up to pressures of ^ 160 GPa. We introduce a 
new set of BM potential parametrizations obtained from 
fits to ab initio energies and atomic forces calculated with 
the van der Waals corrected density functional theory 
method due to Grimme (DFT-D2)ji^ We show that an 
overall improved description of the energy, elastic and 
structural properties of solid helium can be achieved with 
some of the introduced BM three-body interatomic po¬ 
tentials, when used in combination with pairwise poten¬ 
tials in quantum Monte Garlo simulations. Our work also 
brings new insight into the physics of quantum crystals 
at high pressure. For instance, we show that previously 
reported small variations of the kinetic energy, Ek, and 
Lindeman ratio, 7 , in solid helium under pressure^ are 
likely to be artifacts deriving from the use of incomplete 
atomic interaction models. Moreover, we quantify the 
role of quantum nuclear effects on the estimation of the 
shear modulus, Caa, and conclude that they become sec¬ 
ondary when pressure is raised. Finally, at P ~ 25 GPa 
we find that the agreement between our Caa results and 
experiments starts to worsen. Therefore, we argue that 
higher order many-body terms in the expansion of the 
atomic interactions probably are necessary in order to 
describe elasticity in dense solid helium more accurately. 

The organization of this article is as follows. In the 
next section, we outline the employed computational 
methods and provide the technical details in our calcula¬ 
tions. In Sec. imi we explain the fitting strategy that we 
have followed to obtain the new set of three-body inter¬ 
action models. Next, we present our results on the equa¬ 
tion of state, kinetic energy, and structural and elastic 
properties in solid helium, together with some discussion. 
Finally, we summarize our main findings in Sec. [V] 


II. COMPUTATIONAL METHODS 

We used the density functional theory method includ¬ 
ing van der Waals corrections due to Grimme,— to com¬ 
pute the interactions and forces between helium atoms 
in the hexagonal close package (hep) crystal structure, 
from equilibrium up to a pressure of ~ 160 GPa. (De¬ 
tails of our ab initio DFT-D2 calculations can be found in 
elsewhere?^ hence we highlight here only the main tech¬ 
nical features.) Subsequently, we found a series of three- 
body interaction models that, when used in combination 
with the pairwise Aziz-II potential^ (hereafter denoted 
as V 2 ), reproduced very closely the obtained DFT-D2 
results. The details of our fitting strategy are compre¬ 
hensively explained in Sec. IIIII Next, we performed dif¬ 
fusion Monte Garlo (DMG) simulations in which the new 
three-body interaction models were used to calculate the 
energy, structural, and elastic properties of solid helium 
under pressure. In this section, we explain the specific 
implementation of the DFT-D2 and DMG methods in 
our work. 


A. Density functional theory 


We chose the generalized gradient approximation to 
density functional theory proposed by Perdew, Burke, 
and Ernzerhof (GGA-PBE),— as is implemented in the 
VASP package.— Van der Waals interactions were taken 
into account by adding an attractive energy term to 
the exchange-correlation energy of the form Adisp = 
~ Si j C'e /'<’% (where indexes i and j label different par¬ 
ticles, Cq is a constant, and a damping factor is intro¬ 
duced at short distances to avoid divergences)— *21 The 
projector-augmented-wave technique^liH was employed 
to represent the core electrons since this approach has 
been shown to provide very accurate total energies and 
is computationally very efficient— *21 The electronic wave 
functions were represented in a plane-wave basis trun¬ 
cated at 500 eV, and for integrations within the first Bril- 
louin zone (BZ) we employed dense F-centered fc-point 
grids of 14 X 14 X 14. By using these parameters we ob¬ 
tained interaction energies that were converged to within 
5 K per atom. Geometry relaxations were performed by 
using a conjugate-gradient algorithm that kept the vol¬ 
ume of the unit cell fixed and permitted variations of its 
shape. The imposed tolerance on the atomic forces was 
0.005 eV-A“^. With such a DFT-D2 setup we calculated 
the total energy and shear modulus in solid ^He in the 
volume interval 3 < V < 16 A^/atom. 

Additionally, we computed the vibrational phonon 
spectrum in solid ^He at eight different volumes by means 
of the “direct approach”. In the direct approach the 
force-constant matrix is directly calculated in real-space 
by considering the proportionality between atomic dis¬ 
placements and forces when the former are sufficiently 
small.— In this case, large supercells have to be simu¬ 
lated in order to guarantee that the elements of the force- 
constant matrix have all fallen off to negligible values at 
their boundaries, a condition that follows from the use of 
periodic boundary conditions— Once the force-constant 
matrix is obtained, we Fourier-transform it to obtain the 
phonon spectrum at any g-point. The quantities with 
respect to which our DFT-D2 phonon calculations need 
to be converged are the size of the supercell and atomic 
displacements, and the numerical accuracy in the atomic 
forces. The following settings were found to fulfill our 
convergence requirement of correct zero-point energy cor¬ 
rections to within 5 K/ atom)^i2& 4x4x3 supercells (that 
is, 48 repetitions of the hep unit cell containing a total 
of 96 atoms), and atomic displacements of 0.02 A . Re¬ 
garding the calculation of the atomic forces with VASP, 
we found that the density of fc-points had to be increased 
slightly with respect to the value used in the energy cal¬ 
culations (i.e., from 14 x 14 x 14 to 16 x 16 x 16) and that 
computation of the non-local parts of the pseudopoten¬ 
tial contributions needed to be performed in reciprocal, 
rather than real, space. 
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B. Diffusion Monte Carlo 


In our DMC simulations, we used a guiding wave func¬ 
tion, dfsNj, that accounts simultaneously for the atomic 
periodicity and Bose-Einstein quantum symmetry in ^He 
crystals. This model wave function is expressed as^ 

N N / N \ 

'l'SNj(i’l, ■ • ■ ,i"Ar) = 

i<j ,7=1 \i=l / 


where indexes {i,j} and J run over particles and per¬ 
fect lattice positions, respectively. In previous works 
we have shown that provides an excellent descrip¬ 

tion of the ground-state properties of bulk hep "*^He 
and other similar quantum systemsi^^— The correla¬ 
tion factors in Eq. were expressed in the McMil¬ 
lan, /(r) = exp [—1/2 (b/r)^] , and Gaussian, g{r) = 
exp [—1/2 (ar^)], forms. Parameters a and b were op¬ 
timized at each density point by using the variational 
Monte Carlo (VMC) method. For instance, at p = 
0.06 we obtained b = 2.94 A and a = 3.21 A“^ , and 
at p = 0.33 A-3 , 6 = 1.84 A and a = 29.08 A-^ . We 
note that our choice of the guiding function was moti¬ 
vated by an interest in studying the possible effects of 
quantum atomic exchanges on the energetic and elas¬ 
tic properties of dense helium. However, we realised 
by direct comparison to the results obtained with non- 
symmetric wave function models in analogous DMC sim¬ 
ulations?^ that such effects can be totally neglected in 
practice. 

The technical parameters in our calculations were set 
to ensure convergence of the total energy per particle 
to less than 5 K. The value of the mean population of 
walkers was 10 ^ and the length of the imaginary time- 
step (At) 10“"^ K~^ . We used simulation cells contain¬ 
ing 180 atoms. Numerical bias stemming from the finite 
size of the simulation box were minimised by following 
the variational correction approach explained in works y 
and @. Statistics were accumulated over 10® DMC steps 
performed after system equilibration, and the approxima¬ 
tion used for the short-time Green’s function, was 

accurate to second order in Ati^i^ The computational 
strategy that we followed to calculate the shear modulus 
C 44 was the same than in Refs, (ssl - l^ . 


III. FITTING STRATEGY AND THREE-BODY 
POTENTIAL MODELS 


Our three-body potential matching algorithm^^r.^ is 
based on a least square fit to the DFT-D2 reference data, 
that consists of total energies and atomic forces. The 
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FIG. 1. (Top) Energy differences between the DFT-D2 
method and V 3 potentials calculated on a reference set of 16 
configurations (see text). Details are magnified in the inset. 
(Bottom) Results of our fit obtained in the case of the atomic 
forces. AF stands for the difference in the atomic forces be¬ 
tween the DFT-D2 method and many-body potentials, SF for 
the variance of the atomic forces computed with the DFT-D2 
method, and (• • •) for the average performed over particles 
and Cartesian components. 


objective function to be minimized is given by 
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where A" = 16 is the number of reference configurations, 
n = 96 the number of particles on each configuration, and 
loe and luf & weight assigned to the energy, E, and force, 
F, contributions to respectively. With this definition 
of the objective function we ensure that despite differ¬ 
ent magnitudes are expressed in different units all them 
are normalized and contribute equally to x^- Subscripts 
“DFT” and “FF” refer to the DFT-D2 and classical po- 
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FIG. 2. Phonon spectrum obtained with the DFT-D2 method 
(dashed lines) and the V 2 + Vi 3 (BM-F) interaction model (solid 
lines), which was determined considering only the atomic 
forces in the corresponding fit (see text). 


tential results, respectively. 

The set of reference configurations in our fit com¬ 
prised the 16 structures used in the calculation of the 
^He vibrational phonon spectra in the interval 3 < < 

16 A^/atom by means of the “direct approach” (see 
Sec. lII Such atomic arrangements were generated 

by taking the relaxed hep lattice supercells {PG^/mme, 
space group 194) at 8 different volumes and displacing 
one of the atoms sitting in an inequivalent d Wyckoff 
position a distance of 0.02 A first along the ix — 
direction (where x, y, z represent the normalised Carte¬ 
sian vectors), and then along z (that is, we created two 
different atomic configurations at each volume). The rea¬ 
son for our choice was that we wanted to reproduce si¬ 
multaneously the energy and elastic properties in highly 
compressed solid "*^He. In fact, the atomic forces are de¬ 
fined as minus the first derivative of the total energy 
with respect to the atomic positions, whereas the elastic 
constants involve the second derivative of the total en¬ 
ergy with respect to strain deformations. In spite of this 
apparent disconnection, atomic forces and elastic con¬ 
stants are indirectly related by the corresponding spec¬ 
trum of vibrational phonon frequencies. Namely, on one 
side, phonons can be calculated from the variation of the 
atomic forces upon the displacement of atoms away from 
their equilibrium positions, and, on the other side, elas¬ 
tic constants can be estimated from the slope of specific 
acoustic branches in the vicinity of the T point in recip¬ 
rocal space (that is, in the q ^ 0 limit). Therefore, even 
though we did not explicitly consider second derivatives 
in our definition of the objective function we expected 
to achieve an acceptable description of elasticity in solid 
helium. We shall come back to this point later on this 
section. 

The classical potential adopted in this study, denoted 
as “FF” in Eq. ^ is given by Upot = V 2 + V 3 , where V 2 


represents the pairwise Aziz-II interaction model^ and V 3 
the three-body Bruch-McGee (BM) potential given byi^ 


V3 {rij,rik,rjk) = 


1 ^3 - ^exp (-a[ry -b nk + Tjk]) 

' ij ik jk 


X (1 -b 3 cos (j)i cos 4>j cos (j)k) 


(3) 


where =| |, and (j)i, (j)j^ and (pk, are the interior 

angles of the triangle formed by the atoms labelled i, j, 
and k. V 3 is an attractive potential term representing 
triple dipole and three-body exchange interactions. Pa¬ 
rameters V, A, and a were varied during the minimiza¬ 
tion of the objective function (see Eq. [2). Eor this, 
we used a quadratic polynomial interpolation line-search 
with the directions found using the Broyden-Eletcher- 
Goldfarb-Shanno (BFGS) formula4i The gradient of the 
objective function was calculated analytically since other¬ 
wise numerical bias developed that impeded convergence. 
Actually, the typical size of the involved atomic forces is 
very small, of the order of 0.01 — 0.1 eV/A , hence they 
needed to be calculated very precisely. The minimiza¬ 
tions were stopped when all the gradients of the objective 
function in absolute value were smaller than 10“^. Typi¬ 
cally, this was achieved within ~ 100 minimization loops 
when starting from a reasonable initial guess of the A, 
and a parameters (e.g ., the original values proposed by 
Bruch and McGee jl2j|). 

Table I shows the values of the parameters obtained in 
our V 3 fits, in which we considered three different pos¬ 
sibilities based on the choice of the relative energy and 
forces weights: ( 1 ) = 1 and wf = 0 , hereafter de¬ 

noted as V 3 (BM-E), ( 2 ) wf = 0 and wf = 1 , V 3 (BM-E), 
and (3) wf = 0.5 and wf = 0.5, V 3 (BM-EE). Our re¬ 
sults differ appreciably from the original values proposed 
by Brunch and McGee [which hereafter are denoted as 
V 3 (BM)]. For instance, v becomes negative when the en¬ 
ergies are taken into account in the fit, and A and a 
systematically turn out to be larger. 

In Figure [U we illustrate the quality of our fits by plot¬ 
ting the energies and forces calculated on each reference 
configuration. For comparison purposes, we also enclose 
the results obtained with the original V 3 (BM) potential 
(i.e., with function C/pot = V2 -b V3). For the sake of sim¬ 
plifying the notation, we only indicate the three-body 
part in the corresponding many-body potential. This 
convention will be adopted throughout the text if not 
stated otherwise. As is appreciated in the figure, E 3 (BM- 
E) reproduces the DFT-D2 energies more closely than 
any other model (as expected) whereas V 3 (BM) provides 
the worst description. The energies obtained with the 
V 3 (BM-EE) potential can be regarded also as fairly good. 
As for the atomic forces, V 3 (BM-F) produces the best re¬ 
sults, as expected, and E 3 (BM), again, turns out to be 
the less reliable. In this latter case, the forces obtained 
with the V 3 (BM-EE) and, surprisingly also, V 3 (BM-E) 
potentials are not too distant from the reference DFT- 
D2 data. 

In Eig [21 we plot the vibrational phonon spectra ob¬ 
tained with the DET-D2 method and the V 3 (BM-F) po- 
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F (K-cr®) 

A(K) 

a (a ‘) 

V3(BM) [1^ 

0.3270 

9 676 545.53 

4.9480 

V3(BM-E) 

-0.4910 

14 754 161.38 

5.6128 

V3(BM-F) 

1.4029 

12 863 029.73 

5.8273 

V3(BM-EF) 

-1.1364 

29 189 436.37 

6.0691 


TABLE I. Bruch-McGee three-body potential parameters corresponding to the original model, V 3 (BM), and the new ones 
introduced in the present work. The V 3 (BM-E) set has been obtained by considering exclusively DFT-D2 energies on the fit 
[loe = 1, ujf = 0], the V 3 (BM-F) the atomic forces [loe = 0, uf = 1], and V 3 (BM-EF) a combination of ab initio energies and 
atomic forces [loe ~ 0.5, lof = 0.5] (see text). 



3 4 5 6 7 8 9 10 11 12 13 
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FIG. 3. Zero-temperature equation of state calculated in he¬ 
lium with the DFT-D2 and DMG methods. In the DMC 
case, different pairwise and three-body interaction models 
have been employed. Experimental data from Ref. are 
shown for comparison. Inset- The high-P region in the P{V) 
curves are magnified in order to appreciate better the differ¬ 
ences. 


tential in solid "‘He at the highest analysed pressure 
(probably the most challenging case to be reproduced 
with a potential function, see Fig. [I}. We note that the 
agreement between the two sets of data can be regarded 
as fairly good. The largest differences are found on the 
optical branches, which correspond to the highest vibra¬ 
tional frequency values. The DFT-D2 acoustic phonon 
modes in the vicinity of the F point, however, are rea¬ 
sonably well reproduced by V 3 (BM-F). These outcomes 
demonstrate that, as we suggested above, by considering 
the atomic forces in the definition of in principle one 
can obtain a reasonable description of the elasticity in 
the reference system. 


IV. RESULTS AND DISCUSSION 


A. Equation of state 


We show the results of our calculations on the equa¬ 
tion of state, P{V), of solid helium in Fig. [31 together 
with experimental data from work . The DFT-D2 se¬ 
ries was obtained with the ab initio methods explained 
in Sec. Ill A[ including quantum zero-point energy cor¬ 
rections. The other results were obtained in diffusion 
Monte Carlo (DMC) simulations using the indicated in¬ 
teraction potentials, as explained in Sec. Ill HI and else¬ 
where Labels “V' 2 ” and “V' 2 (BC)” stand respectively 
for the pairwise potential due to Azia^ and a modified 
version of the former that we have recently introduced 
in work The DMC (DFT-D2) calculations were per¬ 
formed at 12 (8) different volumes spanned in the interval 
3 < U < 16 A^/atom. In each case, the resulting total 
energies were fitted to a third order Birch-Murnaghan 
equation of the form^i^ 


E(y) -Eo = -VoBoX 


(1 + 2x) 




-^(1 + x) 




(4) 


where Bq = is the value of the bulk modulus at 

the equilibrium volume Vb, x = | ^4 — with Hg = 

(dBo/dP), and all the derivatives are calculated at zero 
pressure. For reproducibility purposes, we enclose the Vb, 
Bq, and Bq parameters obtained in all our fits in Table II. 

Very good agreement is obtained between our DFT-D2 
results and experiments. This outcome justifies in part 
our choice of the DFT-D2 results as reference data in 
modeling of the many-body interactions. Likewise, the 
P{V) curves obtained with the V 2 (BC), V 3 (BM-E), and 
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Vo (A3) 

Bo (eV/A3) 

Bo 

DFT - D2 

12.23 

0.0398 

3.9648 

F2(BC) 

15.68 

0.0166 

4.1144 

V 2 

16.61 

0.0115 

4.8829 

V2-bF3(BM) 

15.68 

0.0181 

3.6722 

1/2 + V'3(BM-E) 

15.84 

0.0165 

4.1854 

1/2 -b V'3(BM-F) 

16.58 

0.0130 

4.6709 

F 2 -b V'3(BM-EF) 

15.85 

0.0158 

4.2463 


TABLE II. Parameters corresponding to the fit of our equation of state results to Birch-Murnaghan functions, see Eq. o , as 
obtained with different computational approaches. In the DMC case, different pairwise and three-body potentials have been 
considered for the description of the interatomic forces. 



FIG. 4. Atomic kinetic energy calculated in “^He with the 
DFT-D2 and DMC methods and expressed as a function of 
pressure. In the DMC case, different pairwise and three-body 
interaction models have been considered for the description 
of the atomic interactions. 


V 3 (BM-EF) potentials are also very close to the obser¬ 
vations. We notice that the V 2 (BC) model introduced 
in Ref. @ was constructed to reproduce the equation of 
state calculated with the DFT-D2 method and that the 
good agreement displayed in Fig. [3] is not a new result. 
Contrarily, the V 2 , 1/3 (BM), and V 3 (BM-F) potentials 
provide a poor description of the variation of the volume 
under pressure. In particular, we find that the V 3 (BM) 
potential systematically underestimates V at pressures 
equal or larger than 20 GPa, in accordance with previous 
results reported by other authors Meanwhile, the V 2 
and Vi 3 (BM-F) interaction models significantly overesti¬ 


mate the same quantity at pressures also close to or larger 
than 20 GPa. In this latter case, we notice a surprising 
resemblance between the two calculated P(V) curves. 

The main conclusion emerging from this part of our 
study is that the new V 3 (BM-F) and V 3 (BM-EF) three- 
body potentials reproduce very accurately the equation 
of state of solid helium up to a pressure of ^ 60 GPa (and 
possibly beyond). To the best of our knowledge, such a 
good agreement between theory and experiments has not 
been reported before for any known V 3 potential in solid 
"‘He (see work [H). 


B. Kinetic energy 

Our kinetic energy, results are shown in Fig. 2] 
In our DFT-D2 calculations, the kinetic energy was esti¬ 
mated within the quasiharmonic approximation through 
the expression 

= ’ ( 5 ) 

^ gs 

where Wqs are the vibrational phonon frequencies in the 
crystal calculated at wave vector q and phonon branch s, 
which depend on the volume, and Nq the total number of 
wave vectors used for integration within the first Brillouin 
zone (see Sec. Ill Al and works @ and [15). usually 
is referred to as the “zero-point energy” (ZPE) and in 
many computational studies turns out to be crucial for 
predicting accurate solid-solid phase transitions 
Regarding our DMC calculations, we computed first the 
exact potential energy, Ep, by means of the pure estima¬ 
tor techniqu o^^d^ and subsequently obtained the exact 
kinetic energy by subtracting Ep to the corresponding 
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total energy. In all the cases, spline interpolations were 
applied to the calculated data points in order to obtain 
smooth P-dependent energy curves (lines in Fig.|4]). 

As is appreciated in the figure, the DFT-D2 results dif¬ 
fer enormously from the rest of Pk series obtained with 
pairwise and three-body potentials in our DMC simu¬ 
lations. At the highest analysed pressure, for instance, 
the DFT-D2 kinetic energy is a factor of two larger than 
the obtained DMC value. Given the lack of experimen¬ 
tal data in the thermodynamic regime of interest, we 
can not rigorously conclude which type of calculation is 
providing the most realistic description. Nevertheless, 
we think that the DFT-D2 results are overestimating 
i?k severely because they have been obtained using the 
quasiharmonic approximation. In fact, it has been al¬ 
ready demonstrated that the quasiharmonic approxima¬ 
tion is not appropriate for studying crystals that behave 
much more classically than solid helium like, for instance, 
molecular hydrogen,—ammonia,—!^ and some alkali 
metals.— It is worth noticing here that although the 
quasiharmonic DFT-D2 approach can produce equations 
of state that are in very good agreement with experi¬ 
ments (as it has been shown in Sec. IIV All , the accompa¬ 
nying ZPE corrections have a lot of margin for error since 
at high P these are always several orders of magnitude 
smaller than the energy of the perfect crystal lattice. We 
shall comment again on this point in the next paragraph. 

It is interesting to analyse the differences found be¬ 
tween the (full quantum) DMC results obtained with 
different pairwise and three-body potential models. The 
V 2 (BC) curve shows a plateau around 550 K at pres¬ 
sures equal and beyond ~ 80 GPa. In a recent work;^ we 
identified such an infinitesimal variation in the kinetic 
energy with the presence of extreme quantum nuclear ef¬ 
fects. However, calculations performed with the new set 
of three-body potentials introduced in this work bring 
new light into our previous interpretation of the V 2 (BC) 
results. As is observed in Fig. SI the Vi 3 (BM-E), Vi 3 (BM- 
F), and V 3 (BM-EF) curves consistently display a small 
but steady increase in the kinetic energy under compres¬ 
sion. At pressures below ^ 15 GPa the pairwise and 
three-body interaction models roughly provide equiva¬ 
lent Ek results however at P = 160 GPa the differences 
between them are as a large of ~ 300 K, with the V 3 
potentials providing always the largest values. Several 
conclusions can be drawn from these results. First, al¬ 
though attenuated pairwise potentials based on exponen¬ 
tial prefactors^ can fairly reproduce experimental P{V) 
data^ii they are likely to introduce unwanted bias on the 
calculation of the kinetic energy. And second, the large 
Ek discrepancies observed between the DFT-D2 and V 3 
results do not seem to be originated by the absence of 
four-, five- and so on many-body interactions in the DMC 
calculations. Actually, by comparing the energy curves 
obtained in the V 2 and V 2 + V 3 cases one realizes that 
the effect of considering three-body interactions on Ek is 
rather small [only in the Vi 3 (BM) case those effects are 
not negligible, although certainly minor]. Therefore, it is 
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FIG. 5. (Top) Atomic density profile around the perfect 
lattice positions calculated with the DMC method consid¬ 
ering different pairwise and three-body interaction models 
(V = 3.0 A^/atom). Solid lines correspond to Gaussian curves 
fitted to the results. The corresponding tails are magnified in 
the inset in order to better appreciate the differences. (Bot¬ 
tom) Lindeman ratio calculated in solid ^He with the DFT-D2 
and DMC methods, expressed as a function of pressure. 

reasonable to expect similar trends when eventually one 
would add higher order many-body terms in the descrip¬ 
tion of the atomic interactions. In regard to this last 
point, we notice that one of the main conclusions pre¬ 
sented in work Q, namely that the quasiharmonic DFT 
approach exceedingly overestimates Ak in dense ^He, ap¬ 
pears to be valid. 


C. Structural properties 

An analysis of the atomic structure in solid ^He at 
high pressure will allow us to understand better the ori¬ 
gins of the discrepancies found so far between the V 2 (BG) 
and V 3 potentials. Figure [5] shows the atomic density 
profiles, ^(r), and Lindeman ratio, 7 , calculated using 
the DMG method and several atomic interaction models. 
The ti{r) results (see top panel) were obtained at volume 





















V = 3.0 A^/atom and subsequently were fitted to Gaus¬ 
sian functions (solid lines in the figure). As is observed 
there, the 1^2(60) curve is noticeably broader than all the 
others, and its value at the origin is about 50 % of that 
calculated with the V3(BM) potential. Meanwhile, the 
V3(BM-E) and V3(BM-EF) profiles are practically indis¬ 
tinguishable and slightly higher near zero than the one 
obtained in the Vi3(BM) case. Clearly, the V2(BC) poten¬ 
tial produces a much larger atomic delocalization than 
the rest of interaction models, which is consistent with 
the kinetic energy results explained in the previous sec¬ 
tion. 

As for the Lindeman ratio 7 (see bottom panel in 
Fig.[S]), we have estimated the corresponding dependence 
on pressure for each analysed potential. In the DFT-D 2 
case, 7 was computed within the quasiharmonic approxi¬ 
mation using the formula see Eq. ([ 5 ]) and 

works [31I and IH^]. The results obtained in the E2(BC) 
case are already known: a plateau around 0.10 appears at 
pressures larger than ^ 80 GPa^ However, all the other 
interaction models, including V2 and V3(BM), provide 
much smaller values of 7 at similar conditions. Moreover, 
the computed Lindeman ratio curves get depleted when 
compression is raised [with the exception of Vi3(BM), in 
which 7 saturates around 0.08 at pressures larger than 
~ 50 GPa]. This latter trend is also observed in the 
DFT-D 2 series, which systematically lies below the DMC 
predictions. 

The results presented in this section show that the 
V2(BG) potential produces an unusually large delocal¬ 
ization of the atoms, which is at odds with the trends 
realised in the rest of cases. Such a huge particle disper¬ 
sion effect is the responsible for the flat kinetic energy 
curve appearing in Fig. 01 which is likely to be an ar¬ 
tifact deriving from the use of exponential attenuation 
factors at short distances. 


D. Elastic properties 

In Figs. El and [71 we show the bulk and shear mod¬ 
uli!, B and C44 respectively, calculated in solid helium 
under pressure. The bulk modulus was directly obtained 
from the Birch-Murnaghan fits explained in Sec. IIV Al 
and in the C44 case spline interpolations were applied 
to the calculated data points in order to obtain smooth 
P-dependent curves. 

Goncerning the analysis of our B(V) results, this is 
very much similar to the conclusions presented for the 
equation of state in Sec. IIV Al Essentially, the DFT-D 2 , 
P2(BG), P3(BM-E), and V3(BM-EF) curves are in good 
agreement with experiments whereas the P2, V3(BM-F), 
and P3(BM) curves are not. In this latter case, both V2 
and V3(BM-F) series are very similar and significantly 
overestimate the bulk modulus at small volumes. Like¬ 
wise, the P3(BM) potential provides unrealistically small 
values oi B{y) at large densities. 

Let us now comment on the (744 (P) results shown in 



V (Adatom) 

FIG. 6. Calculated bulk modulus in “^He using the DFT-D2 
and DMC methods and expressed as a function of volume. 
In the DMC case, different pairwise and three-body interac¬ 
tion models have been considered. Experimental data from 
work 0 are shown for comparison. Inset: The high-P re¬ 
gion in the B{P) curves are magnified in order to appreciate 
better the differences. 



V (A^/atom) 


FIG. 7. Calculated shear modulus in “^He using the DFT-D2 
method and several force fields considering the atoms immo¬ 
bile in the perfect lattice positions. Experimental data are 
from work [ill]. 


Fig. [T] All the values have been obtained considering 
the atoms fixed on their perfect lattice positions, that is, 
totally neglecting likely quantum nuclear effects (hence 
the employed subscript). This is done for the sake of 
comparison since it is technically difficult to account for 
quantun nuclear effects in the DFT-D 2 calculations in an 
exact manner, that is, to go beyond the quasiharmonic 
approximation. Nevertheless, later on this section we will 
show that according to our DMG simulations quantum 
nuclear effects become secondary on C44 at high pressure. 
As is observed in the figure, the DFT-D 2 curve is in over¬ 
all good agreement with the ambient temperature mea¬ 
surements performed by Zha and collaboratorsAi Again, 
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P{V) 

B{V) 

G44(F) 

Ek/7 

General performance 

R 2 [3] 

X 

X 

V/x 

VV) 

Not satisfactory 

F2(BC) 

y 

V 

X 

X 

Not satisfactory 

F2 + F3(BM) [1^ 

X 

X 

Wx 

VV) 

Not satisfactory 

F 2 + V'3(BM-E) 

V 

V 

V/x 

VV) 

Overall good 

F2 + V'3(BM-F) 

X 

X 

VI X 

VV) 

Not satisfactory 

Fa + F3(BM-EF) 


V 

VI X 

VV) 

Overall good 


TABLE III. Summary of the performance of the pairwise and three-body atomic interaction models analysed in this work in 
describing the energy, structural, and elastic properties of solid ^He at high pressure. Symbol V" (^) indicates correct (incorrect) 
description of the considered quantity, whereas V^/x means quantitatively correct up to a certain pressure. Question mark “?” 
denotes a certain hesitation due to lack of experimental data in the high pressure regime of interest. 


these findings justify our choice of the benchmark data 
for the modeling of three-body interactions. Regarding 
the performance of the original and new three-body BM 
potentials, we find that in general they reproduce quite 
satisfactorily the experimental data obtained at volumes 
larger than ^ 5.5 A^/atom (i.e., P < 25 GPa). This is 
especially true in the V 3 (BM-F) case where, as expected 
(see Sec. nn), the calculated shear modulii follow closely 
those obtained with the DFT-D2 method. However, at 
volumes smaller than ~ 5.5 A^/atom (i.e., P >25 GPa) 
we find that the differences between the BM curves (in¬ 
cluding the V 3 BM-F case), on one side, and the DFT-D2 
results and experiments, on the other, become increas¬ 
ingly larger. We recall that the V 3 (BM-E) and P 3 (BM- 
EF) potentials provide a very good description of the 
equation of state and bulk modulus, whereas the V( 3 (BM- 
F) potential does not. This appreciation let us to con¬ 
clude that is very difficult to provide simultaneously a 
good account of the energy and elastic properties in solid 
helium by using an effective three-body approach. Higher 
order many-body contributions in the description of the 
atomic interactions probably are necessary in order to 
attain an overall correct description of solid helium at 
high pressure. As for the pairwise potentials, V 2 performs 
very similarly to the V 3 (BM-F) model, as we have also 
noted in the total energy (see Sec. II V A]) and bulk modu¬ 
lus cases. The V 2 (BG) model, however, remarkably fails 
in reproducing the variation of the shear modulus under 
pressure. Moreover, it predicts the occurrence of unreal¬ 
istic mechanical instabilities (i.e., dC^i/dV ~ at 

small volumes. Therefore, the use of the P 2 (BG) poten¬ 
tial is strongly not recommended for the simulation of 
solid helium at high pressure. 

In order to quantify the importance of quantum nu¬ 
clear effects on the calculation of the shear modulus, we 
carried out additional quantum DMC calculations (see 


Sec. Ill Bl and works for details). To our surprise, 

we found that the quantum and classical shear modulii 
results are very similar. For instance, in the V 3 (BM-F) 
case the Cl^assicai _ ^quantum (where subscript 

“quantum” means calculated with the DMG method) 
amounts only to 2 GPa at P 50 GPa. Similar re¬ 
sults were obtained also in the rest of V 2 and V 3 cases. 
We note that the sign of the differences is always posi¬ 
tive, thus the inclusion of quantum nuclear effects tends 
to lower the classical C 44 values, although in a small frac¬ 
tion (i.e., ~ 5 %). This last finding appears to be con¬ 
sistent with conclusions presented in a recent quantum 
Monte Garlo study by Borda et al.^ in which the ideal 
shear strength on the basal plane of hep "^He was found 
to behave analogously than in classical solids. 

V. CONCLUSIONS 

In Table HI we summarise the performance of the anal¬ 
ysed pairwise and three-body potentials in describing the 
energy, elastic and structural properties of solid ^He at 
high pressure. A number of tips can be drawn from our 
results. First of all, the use of pairwise potentials in gen¬ 
eral is not recommended. These either fail to reproduce 
the equation of state and bulk modulus, i.e., V 2 , or the 
kinetic energy, and structural and elastic features, i.e., 
V 2 (BC), in highly compressed quantum crystals. In this 
context, we urge to employ more versatile many-body 
interaction models. This is the case, for instance, of the 
new three-body BM potentials introduced in this work, 
which represent an improvement with respect to previ¬ 
ously reported similar models. Overall, we recommend 
to consider the V 3 (BM-E) and V 3 (BM-EF) parametriza- 
tions in prospective simulation studies because they pro¬ 
vide the most satisfactory general description of dense 
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solid ^He. Indeed, those interaction models can be safely 
employed, for instance, in atomistic high-P high-T sim¬ 
ulations (either classical or quantum), which are of rele¬ 
vance to planetary sciences. Nevertheless, we must note 
that it remains a challenge to attain a precise description 
of elasticity at high pressure by using effective three-body 
potentials, thus in this latter case consideration of higher 
order many-body terms appears to be necessary. 

Importantly, we have shown that the addition of three- 
body forces corrects for the artificially large atomic delo¬ 
calization found with modified pairwise potentials based 
on exponential attenuation factors. Nevertheless, given 
the lack of structural and kinetic energy measurements 
performed at high pressure, we have not been able to 
quantify the accuracy of our 7 and Pk DMC results 


obtained with the V 3 (BM-E) and V 3 (BM-EF) potential 
models. In this regard, advanced computational studies 
in which both the nuclear and electronic degrees of free¬ 
dom in the crystal were to be treated at the quantum 
level are highly desirable. 


ACKNOWLEDGMENTS 

This research was supported under the Australian 
Research Council’s Future Fellowship funding scheme 
(projects number RG134363 and RG151I75), and 
MICINN-Spain (Grants No. MAT2010-18113, GSD2007- 
00041, and FIS2014-56257-G2-1-P). 


* c.cazorla@unsw.edu.au 
^ jordi.boronat@upc.edu 

^ M. H. Kalos, M. A. Lee, and P. A. Whitlock, Pliys. Rev. 
B 24, 115 (1981). 

^ J. Boronat and J. Casulleras, Phys. Rev. B 49, 8920 
(1994). 

3 R. A. Aziz, F. R. W. McCourt, and C. C. K. Wong, Mol. 
Phys. 61, 1487 (1987). 

^ C. Cazorla and J. Boronat, J. Phys.: Condens. Matt. 20, 
015223 (2008). 

® M. Moraldi, J. Low Temp. Phys. 168, 275 (2012). 

® C. Cazorla and J. Boronat, Phys. Rev. B 91, 024103 
(2015). 

^ T. Omiyinka and M. Boninsegni, Phys. Rev. B 88, 024112 
(2013). 

® D. C. Wallace, Thermodynaimcs of Crystals (Wiley, New 
York, 1972). 

® E. Pechenik, I. Kelson, and G. Makov, Phys. Rev. B 78, 
134109 (2008). 

Y. A. Freiman, S. M. Tretyak, A. Grechnev, A. F. Gon¬ 
charov, J. S. Tse, D. Errandonea, H.-K. Mao, and R. J. 
Hemley, Phys. Rev. B 80, 094112 (2009). 

C.-S. Zha, H.-K. Mao, and R. J. Hemley, Phys. Rev. B 70, 
174107 (2004). 

L. W. Bruch and I. J. McGee, J. Chem. Phys. 59, 409 
(1973). 

M. J. Cohen and J. N. Murrel, Chem. Phys. Lett. 260, 371 
(1996). 

P. Loubeyre, Phys. Rev. Lett. 58, 1857 (1986). 

S.-Y. Chang and M. Boninsegni, J. Chem. Phys. 115, 2629 

( 2001 ). 

C. P. Herrero, J. Phys.:Condens. Matt. 18, 3469 (2006). 
C.-L. Tian, F.-S. Liu, F.-Q. Jing, and L.-C. Cai, J. 
Phys.iCondens. Matt. 18, 8103 (2006). 

S. Grimme, J. Comp. Chem. 27, 1787 (2006). 

J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. 
Lett. 77, 3865 (1996). 

G. Kresse and J. Fiirthmuller, Phys. Rev. B 54, 11169 
(1996). 

C. Cazorla, Coord. Chem. Rev. 300, 142 (2015). 

P. E. Blochl, Phys. Rev. B 50, 17953 (1994). 

G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999). 
C. Cazorla, M. J. Gillan, S. Taioli, and D. Alfe, J. Chem. 


Phys. 126, 194502 (2007). 

S. Taioli, M. J. Gillan, C. Cazorla, and D. Alfe, Phys. Rev. 
B 75, 214103 (2007). 

C. Cazorla and J. Iniguez, Phys. Rev. B 88, 214430 (2013). 
S. A. Shevlin, C. Cazorla, and Z. X. Guo, J. Phys. Chem. 
C 116, 13488 (2012). 

C. Cazorla, D. Alfe, and M. J. Gillan, Phys. Rev. Lett. 
101, 049601 (2008). 

®® D. Alfe, Comp. Phys. Commun. 180, 2622 (2009). 

C. Cazorla, G. Astrakharchick, J. Casulleras, and J. 
Boronat, New Journal of Phys. 11, 013047 (2009). 

C. Cazorla and J. Boronat, Phys. Rev. B 77, 024310 
(2008). 

C. Cazorla, G. Astrakharchick, J. Casulleras, and J. 
Boronat, J. Phys.: Condens. Matter 22, 165402 (2010). 

M. C. Gordillo, C. Cazorla, and J. Boronat, Phys. Rev. B 
83, 121406(R) (2011). 

S. A. Chin, Phys. Rev. A 42, 6991 (1990). 

C. Cazorla, Y. Lutsyshyn, and J. Boronat, Phys. Rev. B 
85, 024101 (2012). 

C. Cazorla, Y. Lutsyshyn, and J. Boronat, Phys. Rev. B 
87, 214522 (2013). 

R. Rota, Y. Lutsyshyn, C. Cazorla, and J. Boronat, J. Low 
Temp. Phys. 168, 150 (2012). 

F. Ercolessi and J. B. Adams, Europhys. Lett. 26, 583 
(1994). 

®® J. Sala, E. Guardia, J. Marti', D. Spangberg, and M. Masia, 
J. Chem. Phys. 136, 054103 (2012). 

M. Masia, E. Guardia, and P. Nicolini, Int. J. Quantum. 
Chem. 114, 1036 (2014). 

J. Nocedal and S. J. Wright, Numerical Optimization 
(Springer, Berlin, 2006). 

P. Loubeyre, R. LeToullec, J. P. Pinceaux, H. K. Mao, J. 
Hu, and R. J. Hemley, Phys. Rev. Lett. 71, 2272 (1993). 

F. Birch, J. Geophys. Res. 83, 1257 (1978). 

C. Cazorla, D. Errandonea, and E. Sola, Phys. Rev. B 80, 
064105 (2009). 

R. Barnett, P. Reynolds, and W. A. Lester Jr., J. Comput. 
Phys. 96, 258 (1991). 

J. Casulleras and J. Boronat, Phys. Rev. B 52, 3654 
(1995). 

G. Geneste, M. Torrent, F. Bottin, and P. Loubeyre, Phys. 
Rev. Lett. 109, 155303 (2012). 



11 


J. M. McMahon, M. A. Morales, C. Pierleoni, and D. M. 
Ceperley, Rev. Mod. Pliys. 84, 1607 (2012). 

M. A. Morales, J. M. McMahon, C. Pierleoni, and D. M. 
Ceperley, Phys. Rev. B 87, 184107 (2013). 

F. Datchi, S. Ninet, M. Gauthier, A. M. Saitta, B. Canny, 
and F. Decremps, Phys. Rev. B 73, 174111 (2006). 

C. J. Pickard and R. J. Needs, Nature Materials 7, 775 
(2008). 

I. Errea, B. Rousseau, and A. Bergara, Phys. Rev. Lett. 
106, 165501 (2011). 


Y. Feng, J. Chen, D. Alfe, X-Z. Li, and E. Wang, J. Chem. 
Phys. 142, 064506 (2015). 

D. A. Arms, R. S. Shah, and R. O. Simmons, Phys. Rev. 
B 67, 094303 (2003). 

G. V. Sin’ko and N. A. Smirnov, J. Phys.; Condens. Matter 
14, 6989 (2002). 

G. Grimvall et ai, Rev. Mod. Phys. 84, 945 (2012). 

E. J. L. Borda, W. Cai, and M. de Koning, Phys. Rev. 
Lett. 112, 155303 (2014). 


